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ABSTRACT 

We present three-dimensional magnetohydrodynamic (MHD) simulations of superbubbles, to study 
the importance of MHD effects in the interpretation of images from recent surveys of the Galactic 
plane. These simulations focus ma inly on atmospheres defined by an exponential density distribution 
and the iDicke v fc Lockmaiil (jTHO) density distribution. In each case, the magnetic field is parallel 
to the Galactic plane and we investigate cases with either infinite scale height (constant magnetic 
field) or a constant ratio of gas pressure to magnetic pressure. The three-dimensional structure of 
superbubbles in these simulations is discussed with emphasis on the axial ratio of the cavity as a 
function of magnetic field strength and the age of the bubble. We investigate systematic errors in the 
age of the bubble and scale height of the surrounding medium that may be introduced by modeling the 
data with purely hydrodynamic models. Age estimates derived with symmetric hydrodynamic models 
fitted to an asymmetric magnetized superbubble can differ by up to a factor of four, depending on the 
direction of the line of sight. The scale height of the surrounding medium based on the Kompaneets 
model may be up to 50% lower than the actual scale height. We also present the first ever predictions 
of Faraday rotation by a magnetized superbubble based on three-dimensional MHD simulations. We 
emphasize the importance of MHD effects in the interpretation of observations of superbubbles. 

Subject headings: ISM: bubbles, magnetic fields - methods: numerical - ISM: individual (W4) 



1. INTRODUCTION 

The combined stellar wind and supernova ejecta of 
groups of O and B stars blow large bubbles in the in- 
terstellar medium. The largest of these bubbles, with 
size scales of 100 pc to 1 kpc are commonly referred to 
as superbubbles. The basic structure of a superbubble 
consists of a hot low-density interior, the cavity, sur- 
rounded by a cool shell of swept-up interstellar medium. 
The continuous formation and dissipation of superbub- 
bles is an important factor in the energy balance of the 
interstellar medium, and determines the locations of dif- 
ferent phases of the inter stellar medium on large scales 
(jMcKee fc Ostrikeilll977f ). Compression of the interstel- 
lar medium in the shell may increase cooling and trigger 
the formation of a new generation of stars. Also, the 
ability of large superbubbles to break out of the disk of 
a galaxy and initiate an outflow of chemically enriched 
plasma and ionizing radiation from the disk into the halo 
has a profound influence on the evolution of galaxies. 

Several examples of well-defined superbubbles have 
been identified in the Galaxy (e.g. Heiles 1984; Ma- 
ciejewski et al. 1996; Normandeau et al. 1996; 
Heiles 1998; Ehlerova & Palous 1999; Callaway et al. 
2000; Reynolds et al. 2001; McClure-Griffiths et 
al. 2002, Pidopryhora et al. 2007). These have 
been studied in detail, thanks to their relative proxim- 
ity. Observations with parsec-scale resolution of neu- 
tral and ionized gas reveal important details about 
the interaction between the hot ejecta and the in- 
terstellar medium. New high-resolutio n surveys of 
Gala ct ic atomic hydrogen (HI ) emis s ion (iTa vlor et al.f 
[200l iMcClure-Griffiths et al.l [20051 : iStilet al. 2006) 
have provided unprecedented images with morpholog- 
ical and kinematic informatio n of Galactic superbub- 
bles (jNormandeau et all 119961 : iMcClure-Griffiths et all 



120031) . Physically interesting parameters are usually de- 
rived from the observations by mea ns of analytic mod- 
els that assum e spherical symmetry (ICastor e t al."1975|: 
IWeayer et al.l Il977i) o r axial symmetry ( Kompan eetj 
' 1960 : !Basuet al.lll999n . 

In this paper we investigate the importance of MHD 
effects on physical quantities derived from observed su- 
perbubbles. The effect of the Galactic magnetic field 
is difficult to model because it introduces anisotropy in 
the medium, which requires three-dimensional numerical 
simulations. The first three-dimensional magnetohydro- 
dynamic s i mulat ions of superbubbles were presented by 
iTomisaka (|1998( ). who discussed the importance of the 
magnetic field i n the break-out o f superbubbles from the 
Galactic disk. iTomisakal (|1998f l also described signifi- 
cant departures from spherical and axial symmetry in 
the shape of a magnetized bubble resulting from the in- 
teraction of the expanding superbubble with the Galactic 
magnetic field. The shape and the size of a superbubble 
depend on the strength and the geometry of the Galac- 
tic magnetic field as much as they depend on the density 
distribution of the ambient interstellar medium. The ex- 
panding superbubble in turn redefines the geometry of 
the interstellar medium and the Galactic magnetic field 
i n a volume se v eral h undred parsecs across. 

iKo rpi et ahl (|1999l) and liie Avillez fc Breitschwerdtl 
(2005) performed thee-dimensional MHD simulations 
that include the evolution of a superbubble in 
a supernova-driven turbulent multi-phase interstel- 
lar medium. The super bubbles in these simula- 
tions show significant departures form symmetry be- 
cause of inhomogeneities in the medium in which 
the super bubble expands . iKorpi et al.l (|1999D and 
Ide Avillez fc Breitschwerdtl ()2005D found that a super 
bubble can break out into the halo in such a medium. 
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Although these simulations provide valuable insight in 
the dynamics of superbubbles in a multi-phase interstel- 
lar medium, it is in general difficult to relate these sim- 
ulations to specific observed super bubbles (see however 
Fuchs et al. 2006 for a simulation of the Local Bubble). 

Astrophysical parameters derived from observations, 
such as the age and the energy of a superbubble, or 
the density distribution of the ambient medium have re- 
lied on symmetric analytic models that do not include 
a magnetic field. The departure from axial symmetry 
imposed by the magnetic field introduces systematic er- 
rors that have not been considered before. In this paper 
we present three-dimensional MHD simulations of super- 
bubbles, and we explore the errors introduced by com- 
monly used methods to determine basic parameters from 
observations. In particular, we study the axial ratio of 
the wind-blown cavity for different times and magnetic 
field configurations. The smaller simulation v olume and 
time s pan used in our simulations compared to lTomisakal 
(jTOOS*) are more tailored to the latitude coverage of the 
Galactic plane surveys. We also calculate the first images 
of Faraday rotation by a magnetized superbubble derived 
from our simulations, and emphasize the importance of 
such simulations to make meaningful predictions in this 
area. 

We describe frequently used analytic models §2 and de- 
tail our numerical setup and methods in §3. Numerical 
results and the analysis of magnetic effects on derived 
parameters are presented §4. The specific case of the 
W4 supperbubble is discussed in §5. Faraday rotation 
by magnetized superbubbles is discussed in §6, and con- 
clusions are presented in §7. 

2. ANALYTIC HYDRODYNAMIC MODELS OF BUBBLES 

The exact shape and size of the bubble depends di- 
rectly on the environment and the strength of the source. 
For weaker sources which produce bubbles whose extent 
are much smaller than the scale height of the gas in 
the galaxy, a solution with an app roximately constant 
density profile may be considered. ICastor etal.1 (|l975l ) 
present a solution for a spherical wind-blown bubble ex- 
panding into a non-magnetized medium, in which the 
radius of the bubble varies according to 
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where Rb is the outer shock of the shell, Lg is the mechan- 
ical luminosity of the source, p is the ambient density of 
the undisturbed interstellar medium, and ig is the age of 
the bubble in Myr. The radius of the contact disconti- 
nuity was found to be OMRb bv I Weaver et all (|1977f l. 

For larger bubbles, associated with powerful sources, 
the situation is more complicated. The gravitational po- 
tential of the Galactic disk produces a density gradient 
per pendicular to the p lane of the disk. At early stages, 
the lCastor et al.l ()1975D solution can still be applied, but 
as the radius of the bubble exceeds the scale height of 



the surrounding medium, the density gradient of the sur- 
r ounding medium affe cts the shape of the bubble. 

iKompaneetsI (|1960D (hereafter K60) proposed a solu- 
tion for the evolution of a point explosion in a non- 
magnetized exponential atmosphere that can be applied 
to this situation. This solution is based on the assump- 
tions that the thermal energy is a constant fraction of 
the energy deposited in the initial blast and that the en- 
ergy distribution is uniform throughout the volume of the 
bubble except close to the shock front where the energy 
density can be two to three times th e mean value (see 
also Bisnovatyi-Kogan & Silich 1995). iBasu et all ()1999l ) 
extended the Kompaneets solution to describe continu- 
ous energy injection, which we use in this paper. The 
source is located at z = in an exponential atmosphere 
with scale height H of the form 



piz) = p„e-^/«. 
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The radius of the bubble in cylindrical coordinates 
takes the following form 



R = 2H arccos 
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is a dimensionless transformed time variable, 7 is the 
adiabatic index, Eth is the thermal energy of the bubble, 
and a. is the volume of the bubble. 

The independent variables are the mass density at 
z = (po); the scale height of the atmosphere (H), and 
the mechanical luminosity of the source (is), which de- 
termines Eth- The unit of time in this model is then 
defined by 



, _ fpoH' 



(5) 



For R — 0, Equation ([3]) reduces to two equations for 
the upper and l ower boundary of the blast wave, zi and 
22 respectively (|Basu et al.lll999D 



zi,2 = -2H\n{ It I 
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As y approaches 2, the top of the bubble in the Kompa- 
neets model expands to infinite height in a finite amount 
of time. Physically, this means the shock acceleration 
in the z-direction becomes infinite because of the strong 
density gradient (|Bisnovatvi- Kogan fc Silichlll995D . The 
bottom of the bubble (z2 , Equation[6|) does not penetrate 
downward more than 2_ff ln2 w 1.4_ff, its location at the 
time of blow-out. Since the top of the Kompaneets model 
reaches an infinite height in a finite amount of time, it 
cannot be a valid solution at later times. However, it 
can provide an adequate solution at early times, if the 
initial conditions are consistent with the assumptions of 
the Kompaneets model, i.e. negligible pressure of the 
ambient medium and negligible inertia of the swept up 
medium. 

With the addition of a magnetic field, no three- 
dimensional analytic solutions exist. In order to capture 
the true evolution of the these bubbles numerical simu- 
lations need to be performed. 
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3. SIMULATION SETUP 
3.1. Goals and Limitations 

We aim to include the physics of the Galactic magnetic 
field in the interpretation of observed superbubbles. The 
derivation of physical quantities from the data, such as 
the scale height of the surrounding medium or the age of 
a superbubble is best served by a model that takes into 
account the magnetic field, but not the complications of 
a preprocessed interstellar medium. We do not include 
the effect of Galactic differential rotation or a Coriolis 
force on the superbubble in our simulations. iTomisakal 
(|1998( ) found that the characteristic time scale for shear 
from differential Galactic rotation is ^--^320 Myr while the 
time scale for shear due to rotation of the bubble by the 
Coriolis force is ~50 Myr. These processes operate over 
significantly longer times scales than those considered 
here (up to 20 Myr). 

The evolution of the bubble will also be affected by 
heating and cooling processes. In this paper we dis- 
cuss adiabatic simulations as an approximation of the 
situation where heating balances cooling throughout the 
life of the superbubble, while we use simulations with 
cooling to explore how the evolution changes if cooling 
dominates over heating. Observations of ionized gas as- 
sociated with superbubble shells indicate that a signif- 
icant fraction of the mass in the shell may be photo- 
ionized by the central star cluster and the surrounding 
interstellar radiation field. The ionized mass in the shell 
of the Orion-Eridaiius sup erbubble is 7 x 10^^^400.5 M© 
([Reynolds fc Ogderil Il979t ). The neutral mass of this 
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iu a4Qo 
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bv iBrown et al.( (|1995D . 



(11976), and 2.5 x 10^rf|oo Mq 

These values indicate that 13 to 28% of the mass of 
the Orion-Erida nus shell is ionized by the stars in- 
side the bubble. iPidoprvhora et all (J2007D found equal 
amounts of neutral and ionized gas in the Ophiuchus 
superbubble, that is possibly ionized from the outside 
as well as from the inside. The s hell of the W4 su- 
perbubble (jNormandeau et al]|1996f ) is clearly visible by 
its thermal radio cont inuum emission and Ha emission 
(jDennison et al.lfl997t l. but not in HI. Photo- ionization 
of a substantial fraction of the mass of the shell indicates 
that heating by photoionization is a non-negligible term 
in the energy budget of the gas in the shell. Simulations 
that include cooling without photo-ionization therefore 
neglect a significant heating term in the energy budget 
of the shell. For computational reasons we cannot solve 
radiative transport in the three-dimensional MHD sim- 
ulations. The absence of cooling in our current simula- 
tions corresponds with the approximation of equilibrium 
between heating and cooling in the shell. In Section [7] 
we discuss the effect of cooling on the axial ratios of the 
cavity. 



3.2. Basic Equations 

Numerical simulations of bubbles have bee n per- 
forme d by many resear c hers (see iRobinso n et al.l 

(I2004D. iKomlienovicetal] (Il999f). I^misakal (I1998D . 
Mineshige et all (I1993D . iMac Low et al.l (119891) and 
Tomisaka fc Ikeuchil (|1986l ). among others). ITomisakal 



TABLE 1 

DiMENSIONLESS PARAMETER CONVERSION 



Parameter 


Variable 


Conversion 


Density 


P 


PPo 


Location 


X 


xH 


Velocity 


V 


VCs 


Time 


t 


iH/cs 


Luminosity 


L 


LpoH^cl 


Pressure 


P 


PPocl 


Internal Energy 


e 


epocl 


Magnetic Field 


B 


^P-'I^pT-s 




TABLE 2 




Input Parameters 


Atmosphere 




Source 


PO 




Ls 


/3o 




Vs 


B:(Bi,B2,B3)'' 




Rs 


Density profile'' 


f^on 






^off 




Location: (xi ,a;2 ,X3 ) 



'^ In the simulations presented here the ini- 
tial magnetic field is oriented along the xi 
axis (B2 = S3 = 0), parallel to the Galac- 
tic plane.'' Exp = Exponential according to 
EquationlH DL = Dickey & Lockman (1990) 
according to Equation 1151 both with density 
gradient along the X2 axis 

(|1998f ) was the first to perform three-dimensional MHD 
simulations with radiative cooling assuming symmetry 
with respect to the x=0, y=0 and z=0 planes. In the 
simulations we present here, we have a complete three- 
dimensional bubble evolving in an unperturbed, three- 
dimensional environment. 
Our simulations solve the following equations: 



dp 

dt 



-HV-(pv) = 



dB 

'dt 



dv 

- + (v.V)v 



- V X (vx B) = 



(V X B) X B 



47r 







(7) 
(8) 

(9) 



- + (v.V)e 



+ p(V-v) = (10) 



V-B = 



(11) 



where p is the density, p is the gas pressure, B is the 
magnetic field, v is the velocity, <I> is the gravitational 
potential, and e is the internal energy of the gas. 

For our simulations we work in units of density (po)j 
scale height (if) and sound speed (c^). The conversion 
to dimensionless variables takes the following form 



PPo. 



(12) 



The other parameters are converted in a similar manner 
(see Table [Ij . The unit of time in our simulations is 



H 

to = — 



(13) 
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which is different from the unit of time in the Kompa- 
neets model io,fe in Equation (O, which assumes Cs — 0. 
The MHD equations in dimcnsionless vari ables are 
solve d numerically by the ZEUS-MP code (jNormanl 
|2000( ). The ZEUS-MP code used in this work is an 
adapted and modified version of the original 1.0b ver- 
sion of Norman and Li (see Norman 2000). One of the 
authors involved in this paper (R. Ouyed) has also spent 
a large amount of time to modify the code. Most of 
these modifications w ere similar to those described in 
IVernaleo fc Revnoldsl ([2006) • A detailed description of 
the original ZEUS code can be found in Stone & Nor- 
man (1992a&b) which also includes basic tests of the 
code. Vulnerabilities in the ZEUS family noted by iFalld 
(|2002f ). in particular the issue of rarefaction waves and 
shock errors, have recently been discussed by Hayes et 
al. (2006; and references therein). While these papers 
acknowledge that the code possess limitations as do all 
numerical methods, results from ZEUS-MP were found 
to compare quite favorably with other numerical tech- 
niques. We come to the same concl usion in jjRi.S. 21 w here 



we compare our results to that of iTomisakal ()1998l ) de- 
spite differences in numerical methods used. 

ZEUS solves the MHD equations using the operator 
split method. The equations are solved in two sub- 
steps, called a source step and a transport step (see 
Stone & Norman 1992a&b). Three methods for the 
advection of mass, momentum and internal energy in 
the transport step can be impl emented : the first order 
accurate donor cell method (Godunov '1959"), the sec- 
ond order van Leer method (Van Leer 1977), and the 
third order piecewise parabo lic advection (PPA) method 
(jColella fc WoodwardI [l98^ . After some numerical ex- 
periments we decided to use the van Leer method because 
it offers the best ratio of precision to computational costs. 
The basic equations of the code are written in a covariant 
form which allows for the use of the code in an arbitrary 
orthogonal coordinate system (Cartesian, cylindrical and 
spherical coordinates are predefined). The algorithm 
used to guarantee that V • J5 = is the "HSMOCCT" 
method w hich combines the constrained transport (CT) 
module of lEvans &i HawlevI (|1988D . and improvements 
of the method o f char acteristics (MOC) introduced by 
iHawlev fc Stond (119951 ). In this scheme, if the initial B 
has zero divergence in the discretisation on the staggered 
mesh, then every time step will maintain the initial value 
of the divergence to the accuracy of machine round-off 
error. 

The finite difference method is based on the discretiza- 
tion of each dependent variable over the spatial compu- 
tational domain. Then finite difference approximations 
to the differential equations are solved on this discrete 
mesh. The ZEUS code uses a staggered mesh built up 
of two mutually shifted grids. The a-grid specifies po- 
sitions of the zone boundaries while the b-grid specifies 
the zone centers. Discrete values of all dependent vari- 
ables are stored for each zone. Scalars are stored at the 
zone centers while components of vectors are stored at 
the appropriate zone interfaces (see Figure 1 in Stone & 
Norman 1992a). Boundary conditions are implemented 
as two layers of ghost zones at each boundary of the com- 
putational domain (two layers are required for higher or- 
der interpolation if the PPA method is used) . Values of 
the dependent variables in the ghost zones are given by 



simple, explicit equations that connect these values to 
the values in the adjacent active zones. For our simula- 
tions we use the so-called outflow boundary conditions 
where the values of all variables in the ghost zones are 
set equal to the values in the corresponding active zones. 
The simulations are performed in a three-dimensional 
cartesian box with right handed coordinates xi, X2, and 
X3. As initial conditions, we specify a density distribu- 
tion (atmosphere) and the location of one or more energy 
sources and their luminosities (as described in the follow- 
ing sections) . For a summary of these input parameters 
refer to Table H 

3.3. Atmosphere Setup 

We consider two functional forms for the distribution of 
the ambient gas. The first form is an exponential density 
distribution 

p = exp [-X2] ■ (14) 

The second form is the den sity distribution proposed 
bv lDickev fc LockmanI (|1990D (hereafter DL) from their 
analysis of the vertical distribution of atomic hydrogen 
in the Galaxy. This density distribution, in terms of 
dimcnsionless variables, is 



0.395 
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0.566 



exp 



-03^2! 



(15) 



where ai = 100/90, 02 = 100/225, and 03 = 100/403. 
The coefficients ai express the scale heights of the three 
components of the DL layer in units of 100 pc. Contrary 
to the exponential profile, the DL layer has an equatorial 
plane, p{x2 = 0) = 1, with density decreasing in both 
the positive and negative X2 dir ections. This is t he same 
density distribution adopted bv ITomisakal (|1998( l. 

The dynamical importance of the magnetic field is set 
through the parameter f3a, defined as the ratio between 
the gas pressure and magnetic pressure at 0:2 = 



(^0- -fw 



Snp 

"B2 



(16) 



The magnetic field strength is 



V/3o ; 
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1.67 X 10-24 g cnr 
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where the numeric expression was derived using cj. — 

IP/P- 

We consider two geometries of the magnetic field. The 
first is a constant field, i.e. the vertical scale height of the 
magnetic field is infinite {B{x2) = -Bo), the other takes 
the scale height of the magnetic pressure to be equal to 
the scale height of the gas pressure (/3 = /3o everywhere) . 
We refer to the constant /3 case as equipartition although 
strictly speaking equipartition implies (3=1 everywhere. 
The scale height of the Galactic magnetic field is not well 
known. Our initial conditions cover the possibility of a 
large scale height > 1 kpc corresponding to a magneti- 
cally dominated halo, and a small scale height (equal to 
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Fig. 1. — The oscillation of the source in xi, X2, and X3 results in 
a nearly isotropic source with radius only a few zones. The dashed 
circle represents the idealized spherical source at one time, the solid 
circle represents the idealized source at another time. Horizontal 
and vertical lines represent zone boundaries. The two arrows illus- 
trate velocity vectors for two locations on the a-grid where mass 
outflow from the source occurs when the idealized source is repre- 
sented by the dashed circle. In the simulations, the amplitude of 
the oscillation is set to only half a zone. 

2H for an exponential atmosphere). Hydrostatic equi- 
librium is assured by imposing a gravitational potential 
that balances the gradient of the total pressure which is 
the sum of the gas pressure and the magnetic pressure. 
The maintenance of hydrostatic equilibrium was tested 
in simulations with no source for both atmospheres. The 
vertical velocity in these tests was less than 10~^Cs after 
10 Myr simulated time. This is consistent with numerical 
errors and it has no effect on the results of our simula- 
tions. 

3.4. The Energy Source 

The source of the bubble is defined by its mechanical 
luminosity (Lg) related to the mass loss rate (m^) and 
the outflow velocity (vg) according to 



Ls = ■^rhsv'^g. 



(18) 



The source is assigned a finite radius (Rs), determined 
by the spatial resolution of the simulation. The mass loss 
rate from the source is related to the outflow velocity and 
the radius of the source via 



rhs = psVs^TrRg 



(19) 



Solving for ps and converting to dimensionless vari- 
ables, we find 

Ps = -^ (20) 

Values for Kg, Lg and Vs are specified as initial con- 
ditions. At every time step, ps and Vs are reset to their 
initial values to maintain constant momentum and en- 
ergy input from the source. To simulate a point source, 
the radius Rs should be taken as small as possible. How- 
ever, in the rectangular grid, the outflow from the source 
will be more isotropic if Rs is more than a few times the 
grid size. A complication is that density is evaluated in 
the center of a pixel (b-grid), whereas velocity is evalu- 
ated on the boundary of a pixel (a-grid). This is part 
of a strategy to conserve mass in the simulations, but it 
introduces an asymmetry in the radius of a sphere when 
coordinates are rounded off to integer multiples of the 



TABLE 3 

Physical Parameters for Simulations Performed: 

Source and Atmosphere 



Simulation 


Ls 
(10»^ erg s 


no (z = 0) 
-1) (cm-3) 


B 


/3o 


Constant Atmosphere 


ConstH_zoom 


3 


1 


noB 


cx> 


ConstH 


3 


1 


no B 


CX) 


ConstMHD 


3.2 


0.32 


const 


1.16 


Exponential Atmosphere 


ExpH_zoom 


3 




no B 


00 


ExpH 


3 




no B 


00 


ExpCBa 


3 




const 


10 


ExpCBb 


3 




const 


3 


ExpCBc 


3 




const 


1 


ExpCBd 


3 




const 


0.3 


ExpEBa 


3 




equip 


10 


ExpEBb 


3 




equip 


3 


ExpEBc 


3 




equip 


1 


ExpEBd 


3 




equip 


0.3 




DL 


Atmosphere 






DLH 


3 




no B 


00 


DLCBa 


3 




const 


10 


DLCBb 


3 




const 


3 


DLCBc 


3 




const 


1 


DLCBd 


3 




const 


0.3 


DLEBa 


3 




equip 


10 


DLEBb 


3 




equip 


3 


DLEBc 


3 




equip 


1 


DLEBd 


3 




equip 


0.3 


Tomisaka A 


3 


0.3 


const 


0.7 



grid size. To maintain an isotropic source while avoid- 
ing computationally intensive interpolations that must 
be performed every time step, the location of the source 
is oscillated in three dimensions by half a grid position 
every time step (see Figure [T]) . These excursions of the 
source along the three axes are incoherent, and average 
out most of the anisotropy introduced by the numerical 
grid for small Rs . The isotropy of the source setup was 
tested in hydrodynamic simulations with no density gra- 
dient. These experiments showed that Rg between two 
and three pixels results in good symmetry in the simula- 
tions. A small degree of anisotropy in the outflow of the 
source remains. This acts as a seed for the development 
of instabilities, but is not considered a problem for our 
analysis. 

4. RESULTS 



The simulations were performed on the CAPCA"'^ com- 
puter cluster which consists of 64 2.4 GHz Linux based 
processors connected by a 1 gigabit network. The re- 
s ults of the se simulations were analyzed using JETGET 
(Staff etal. 2004) and KVIS, a p rogram that is part of 
the Karma visuahzation package ()Goochlll996[ ). 

In the simulations that we present here, we exam- 
ined the effects of magnetic field strength and geome- 
try on the evolution of a superbubble in the two atmo- 
spheres defined in Section 13.31 The other parameters 
were fixed to the values in Table [3] to apply our sim- 
ulations to published results on the bubble associated 
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with t he ^ y4 region a s desc ribed by iNormandeau et al.1 
(|1996l ) and lBasu et"al] (|1999l ) (SectionE]). Table[3]gives a 
summary of the simulations and the values of the varied 
physical parameters. Table[4]gives details on the simula- 
tion volume. To test convergence we ran simulations at 
100^, 200'^ and 300'^ zones and established convergence at 
200^ which corresponds to a resolution of 5 pc per voxel 
(see Section |4XT|) . 



4.1. Comparison to Analytical Solutions 
4.1.1. Hydrodynamic Solutions 

In order to test the setup, we performed hydrody- 
namic simulations with a constant density profil e (Con- 
stELzoom & ConstH) to compare with the iCastor et al.l 
(|1975f ) model and a hydrodynamic simulation with an ex- 
ponential density profile (ExpH_zoom) to compare with 
the K60 model. The physical parameters for these simu- 
lations are given in Table [3l Maintaining the same num- 
ber of pixels for a smaller volume (300 pc on a side) 
results in a finer resolution (1.5 pc compared to the 5 
pc resolution of the other simulations listed in Table |4|) . 
This allows us to follow the evolution of the bubble in de- 
tail at early times. We have multiple reasons to consider 
high-resoluti o n sim ulations for the comparison with the 
ICastor et al.l ()1975f ) and K60 models. The first reason 
is to investigate the effect of resolution on our solutions. 
The second reason is to explore the effect of the finite 
source size. In our simulations, the bubble begins at a 
finite radius at t = (the radius of the source) which can 
be made smaller in high resolution simulations. Although 
the ICastor et al.l (|1975( ) model is self-similar, increasing 
the resolution will allow us to investigate the effect of the 
finite source size in the simulations at early times. The 
K60 model is not self-similar, but it is only expected to 
agree with the simulations at early times (Section |2|) . 

Both models assume that the sound speed in the am- 
bient medium is zero by assuming the pressure in the 
undisturbed medium is negligible. For our simulations, 
the sound speed in the undisturbed medium is finite, but 
the assumption of negligible pressure is satisfied as long 
as the pressure inside the bubble is much larger than 
the pressure of the surrounding atmosphere. The limit 
of negligible ambient pressure adopted in the K60 model 
implies to 3> to^ ■ In our simulations, the time scale to 
(Equation lT3| is ^ 10^ years, whereas the timescale to^. of 
the K60 model (Equation JS]) is -- 0.05 x lO"^ years. The 
simulations should therefore resemble the Kompaneets 
model at times much less than 10^ years. 

The size of the bubble was parameterized by the ra- 
dius of the contact discontinuity, which is well defined in 
both the simulations and the Kompaneets model. The 
shell of swept-up interstellar medium is presumed to be 
infinitely thin in the Kompaneets model, so the radius 
of the outer shock is not predicted. For consistency, the 
same was done for the spherical case. The radius of the 
contact discontinuity for the spherical mod el was taken to 
be 86% of the outer shock radius following I Weaver et al.l 
()1977l ). For ConstH_zoom, the contact discontinuity is 
located at 84% of the out er shock radi u s, wh ich is in 
good agreement with the I Weaver et al.l (|1977f ) model. 
Figure [HA compares the time evolution of the radius of 
ConstH_zoom with the spherical model. The simulation 



displays a power law expansion with time i?b ^ i" with 
a = 0.566. This i s only 5.7% smaller than the slope of 
the ICastor et al.l (|1975D model, which has a = 0.6. At 
an age of 1 Myr the radius of ConstH is 9% larger than 
the radius of the model. The filled squares in Figure [21^ 
show the time evolution of ConstH. The radius of the 
bubble in this simulation is at most 10% larger than the 
radius in the higher resolution simulation at the same 
age. This indicates that part of the difference between 
our simulations and the Castor model may be related to 
the finite size of the source. 

Figure |2j3 shows the time evolution of the radius^ of 
a bubble in an exponential atmosphere (ExpH_zoom in 
Tables O and m , looking perpendicular to the Galactic 
plane, for the simulations and the K60 model. Also 
shown is the spherical model from Figure |2j'V as a ref- 
erence line. The bottom of the simulation closely follows 
the K60 model. However the top of the K60 model is 
consistently higher than the simulation by about 16% 
even at early times. This is also shown in Figure [31 As 
time proceeds, the difference between the K60 model and 
the simulations increases sign ificantly at the top of the 
bubble but not at the bottom. iKomljenovic et al\ ()1999l ) 
also compared their 2-dimensional hydrodynamic simu- 
lations to the Kompaneets model with similar results. 
iMac Low et al.l (1989) found consi stency between their 
Komp aneets approximation from iMac Low fc McCravl 
(|l988i ) and hydrodynamic simulations at times as late 



as 6.87 Myr. However, their Kompaneets approximation 
is a solution in an atmosphere with a n equatorial plane 
unlike the Kompaneets model (|Kom panccts 1960) con- 
sidered here. The axial ratio in the Galactic plane (i.e. 
X1X3) is expected to be unity at all times for the hydro- 
dynamic case and also at early times in the MHD case 
(see Section IT^ . The radius of the cross-section through 
the source of the K60 model as a function of time can 
be calculated to a good approximation from the spher- 
ical C astor model, as shown numerically by iBasu et al.l 
(fl999l) . 

4.1.2. MHD Analytical Approximation 



iFerriere et al.l ()199lD derived analytic solutions of a su- 
perbubble in a uniform magnetic field in the limit of high 
expansion velocity, and numerical solutions for the gen- 
eral case of smaller expansion velocity. To compare with 
their solutions, we ran an MHD simulation in an atmo- 
sphere with_const£mt densi ty (ConstMHD, see Tables |3I 
&[4l). IFerriere et al.l (|1991h found that the outer shock 
front remains nearly spherical in a uniform magnetized 
medium, but the cavity is smaller in the direction per- 
pendicular to the magnetic field, compared with the hy- 
drodynamic solution. We find the same in our numerical 
simulation, but we see no dimple in the outer sho ck in 
the direction of the magnetic field. iTomisakal (|1998f ) also 
did not see this dimple. The expansion along the mag- 
netic field lines is nearly the same as in the hy drodynamic 
simulation as found by Ferriere et al.l (|1991h . 

The thickness of the shell perpendicular to the mag- 
netic field in the analytic approximation found by 

^ Here, the radius is defined separately for eacfi direction as 
the distance from the source to the contact discontinuity at the 
top of the bubble and the distance from the source to the contact 
discontinuity at the bottom of the bubble. 
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Box Parameters for Simulations Performed: Volume 


setup" 






Simulation 


a;i X 12 X x-j 


{xiSzX3)min-max {^2 


)min — max 


n^^ X r 


^2 > 


nj;3 




pc X pc X pc 


H 


H 


zones X zones 


X zones 


Constant Atmosphere 


Const ii_zoom 


300 X 300 X 300 


-1.5 to 1.5 


-1 to 2 


200 X 


200 


X200 


ConstH 


10^ X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ConstMHD 


10^ X 103 X 103 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


Exponential Atmosphere 


ExpH_zoom 


300 X 300 X 300 


-1.5 to 1.5 


-1 to 2 


200 X 


200 


X 200 


ExpH 


10» X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpCBa 


10^ X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpCBb 


10^ X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpCBc 


10^ X 103 X 103 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpCBd 


10^ X 103 X 103 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpEBa 


10^ X 103 X 103 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpEBb 


10^ X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpEBc 


10^ X 10^ X 10^ 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


ExpEBd 


103 X 103 X 103 


-5 to 5 


-3 to 7 


200 X 


200 


X 200 


DL Atmosphere 


DLH 


10^ X 10^ X 10^ 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLCBa 


10^ X 10^ X 10^ 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLCBb 


10^ X 10^ X 10^ 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLCBc 


103 X 103 X 103 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLCBd 


103 X 103 X 103 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLEBa 


10^ X 10^ X 10^ 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLEBb 


103 X 10^ X 10^ 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLEBc 


103 X 103 X 103 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


DLEBd 


103 X 103 X 103 


-5 to 5 


-5 to 5 


200 X 


200 


X 200 


Tomisaka A 


10^ X (2 X 103) X IC 


3 to 10 


-10 to 10 


200 X 


400 


X 200 



^ This tabic gives information on the main set of simulations. Additional simulations for convergence tests 
with 100 X 100 X 100 and 300 X 300 x 300 zones (see Figure [To} are not listed. 
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Fig. 2. — (A) Time evolution of the radius of a stellar wind bubble in a medium with constant density and no magnetic field. Diamonds 
show the radius of the bubble in a simulation with 200 X 200 X 200 zones that measures 300 X 300 X 300 pc (i.e. a resolution of 1.5 pc per 
voxel). Squares show the radius of a bubble in a 200 X 200 X 200 simulation that measures 1000 X 1000 X 1000 pc (i.e. a resolution of 5 pc 
per voxel). In both cases the rad i us of the source was 2 pixels. The radii of the simulations are compared with the radius of the contact 
discontinuity in the [Weaver et aT] I I1977I) model (dashed line). Also shown is the radius of the outer shock in the Weaver et al. ( 1977) model 
(solid line). (B) Time evolution of the radius of the bubble in an exponential atmosphere in the positive X2 direction (filled triangles) 
and the negative X2 direction (stars). The correspo nding radii from the K60 model are shown as a dashed curve and a dot-dashed curve 
respectively. Also shown is the [Weaver et ahl I I1977I) constant density model from panel A (solid line). See Appendix A for a discussion on 
the difference between the simulations and the K60 model at early times. 
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Fig. 3. — Early stages of a zoomed-in hydrodynamic simulation with Koinpaneets solution overlaid (white circles) in a plane through 
the centre of the source at (xi,a;2,X3)=(0,0,0). Gray scales represent density on a logarithmic scale from 10~^ to lO^'^ cm~^. At early 
times, the Kompaneets model is nearly spherical, but the center of the sphere has a visible offset in the direction of the density gradient 
(see Appendix A). The geometric center of the Kompaneets solution as defined in Appendix A is shown as a + in each panel. The values 
for y (Equation IJjl are 0.36, 0.55, 0.70, and 0.84 for these panels in increasing time order. The unit of length on the axes is 100 pc. 
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iFerriere et al.l (|1991h . increases linearly with time accord- 
ing to 

■ smt' pc (21) 



AR= 1.2- 



where no is the number density in cm^'^, B^q is the am- 
bient magnetic field strength in /iG, and ig is the time 
in Myr. Our simulation reproduces the linear increase of 
shell thickness perpendicular to the magnetic field with 
time, but the s hell thickness is a factor ~ 2.6 larger. The 
simulations by iTomisakal ([19980 ^-l^o show a very thick 
shell perpendicular to the magnetic field, although he did 
not consider a medium with constant density. The rea- 
son for the difference is that Equation (f2T|) was derived 
in the limit of high expansion velocity, defined as a small 
ratio of magnetic pressure to ram pressure. This ratio 
increases with time as the outer shock decelerates. We 
find that this ratio is of order unity or larger at times 
> 2.5 Myr in our simulation. Simulations with higher 
f3 reproduce the shell thickness in Equation pTj) . while 
also satisfying the assumptions made in the derivation of 
this equation. 

4.2. Exponential Atmosphere 
4.2.1. Model Exp H 

A simulation of a bubble expanding in an exponential 
atmosphere with no magnetic field was done primarily to 
compare with previously published results and to obtain 
the limiting case j3 = oo. At early times, the bubble re- 
mains spherical (Figure [3]) . The axial ratio of the cavity 
in the simulations and the axial ratio of the Kompaneets 
model are indeed found to be very close to unity, even 
though Figure [3] shows that the top of the Kompaneets 
model has proceeded significantly further than the sim- 
ulation even at these early times. In Appendix A we 
show that the Kompaneets model at early times can be 
described to third order in j/ as a spherical cavity that 
rises in the atmosphere. Although the axial ratio remains 
unity to a high degree of accuracy, the geometric center 
of the model, defined by the point midway between the 
top and the bottom (crosses in Figure [3]), is displaced 
significantly from the location of the source. 

At later times, the bubble grows larger than the scale 
height of the medium and accelerates in the vertical di- 
rection. As the bubble expands into the upper atmo- 
sphere, the top of the bubble continues to accelerate 
and a Rayleigh- Taylor instability develops when the top 
reaches ^ 4 scale heights above the source, which occurs 
at an age of ~10 Myr. The bubble breaks out of the 
galactic plane to form a chimney after about 16 to 17 
Myr. 

4.2.2. Models ExpCB(a-d) & ExpEB(a-d) 

The main difference between a constant magnetic field 
and equipartition magnetic field is in the evolution of 
the bubble at large distances from the Galactic plane. A 
constant magnetic field implies a high Alfven speed in 
the Galactic halo whereas constant (3 implies that the 
Alfven speed in the halo is the same as in the disk. This 
difference has significant consequences for the structure 
of the shock as it expands from the disk into the halo. 
Figure [4] shows a bubble in an exponential atmosphere 
with constant magnetic field (/3o = 1) oriented along the 
xi axis, at an age of 10 Myr. 



The cavity is elongated along the magnetic field as the 
plasma can move freely along the field lines but its mo- 
tion perpendicular to the field lines is restricted. How- 
ever, the shape of the outer shock at the level of the 
source in the Galactic plane {xix^ plane) is nearly cir- 
cular because the propagation speed of the outer shock 
is similar in both directions. Along the xi axis, the 
outer shock and the contact discontinuity are close to- 
gether resulting in a thin compressed shell of swept-up 
interstellar medium. Along the x^ axis, perpendicular 
to the initial magnetic field, the distance between the 
outer shock and the contact discontinuity is much larger 
resulting in a thick shell. Looking along the field lines, 
the cavity is more elongated than in the hydrodynamic 
simulation, and the shell of swept-up interstellar medium 
is thicker than in the hydrodynamic case. At the top of 
the cavity a fast magnetosonic wave runs upwards into 
the halo, with very little compression of the halo gas. In 
contrast to the hydrodynamic simulation, we see no evi- 
dence for a Rayleigh- Taylor instability developing at the 
top of the bubble, because the top of the bubble does 
not accelerate. Instead, an instability develops every- 
where at the contact discontinuity as is apparent from 
the scalloped shape of the contact discontinuity in the 
x^X2 plane (looking along the field lines). This instabil- 
ity also occurs in the simulations of Tomisaka (1998) and 
may be similar to the magnet ically enhanced Raylcigh- 
Taylor instability proposed by iGregori et al.l (2000. ) . At 
the time shown in Figure HI the expansion of the cavity 
along the 2:3 axis (perpendicular to the magnetic field) 
in the Galactic plane has stalled as a result of magnetic 
tension. The cavity still expands along the xi and X2 
axes. 

Figure [5] shows simulation ExpEBc also at the age of 
10 Myr. In contrast with Figure [H Figure [5] shows a 
shell of compressed interstellar medium at the top of the 
bubble. Since the Alfven velocity does not increase with 
distance from the Galactic plane, the outer magnetosonic 
shock does not accelerate into the halo but remains close 
to the contact discontinuity. The shell appears approxi- 
mately twice as thick looking along the field lines as look- 
ing perpendicular to the magnetic field. The axial ratio 
of the cavity in the X3X2 plane (looking along the field 
lines) is smaller than in the ExpCBc simulation. This 
indicates the importance of magnetic tension confining 
the bubble, even though the total pressure mimics the 
pressure in the hydrodynamical exponential atmosphere. 
The magnetic field in the upper part of the shell is signif- 
icantly enhanced compared to the magnetic field in the 
undisturbed atmosphere at the same height above the 
plane. Here too, we see an instability everywhere along 
the contact discontinuity in the x^X2 plane. 

4.3. Dickey & Lockman Atmosphere 
4.3.1. Model DLH 

The simulations in the DL atmosphere show some dif- 
ferences compared to the evolution in the exponential 
atmosphere. The equatorial plane introduces symmetry 
with respect to the plane 2:2 = 0. In the absence of a 
magnetic field the bubble is spherical until it expands to 
radius ~ 200 pc. Once the radius of the bubble grows 
beyond ~200 pc, it becomes more elongated in the ver- 
tical direction as it balloons out into the halo (Figure [6]). 
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Fig. 4. — Simulation ExpCBc at an age of 10 Myr. Magnetic field initially oriented along xi axis, [3=1. Panels show slices through the 
cube at the location of the source, (3;i,a;2,a;3)=(0,0,0), in three orthogonal planes. Grayscales show the gas density on a logarithmic scale 
from 10~'' (white) to lO^'^ (black) cm~^. The vector field depicts the projection of magnetic field vectors on the plane of the image. The 
unit of length on the axes is 100 pc. 
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Fig. 5. — Same as Figure l4l but for Simulation ExpEBc at an age of 10 Myr. 
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The top and bottom of the bubble expanding into the 
low-density outer atmosphere, accelerate away from the 
Galactic plane, as in the case of an exponential atmo- 
sphere. The onset of a Ray leigh- Taylor instability is not 
seen in this simulation but it may develop at later times 
outside the simulation volume. 

Between an age of ~ 14 Myr and ~ 18 Myr a small de- 
crease in the radius of the cavity in the equatorial plane 
is observed. A similar stall of the expansion and subse- 
quent contraction in the Galactic plane triggered by the 
rapid vert i cal ex pansion of the bubble was observed by 
lTomisakal(|1998f ). 

4.3.2. Models DLCB(a-d) & DLEB(a-d) 

Figure [7] shows the simulation of a bubble in the DL 
atmosphere with constant magnetic field oriented along 
the xi axis at an age of 10 Myr. As with the exponen- 
tial simulations, this simulation s hows the same features 
as discussed by iTomisakal ()1998f ). The cavity is elon- 
gated along the magnetic field lines. The shape of the 



outer shock at the level of the source in the 0:1X3 plane 
is nearly circular, similar to the exponential case. This 
is to be expected because the density and magnetic field 
strengths are the same at the level of the source in both 
atmospheric profiles. In contrast to the exponential at- 
mosphere, the shape of the cavity in the xia;2 plane is 
fairly circular even though the bubble expands to a size 
much larger than the pressure scale height. The shape of 
the cavity in the 0:3X2 plane is quite elongated and dis- 
plays the same magnetically enhanced Rayleigh- Taylor 
instability as discussed in Section 14.2.21 As with the 
ExpCBc case, a fast magnetosonic wave runs upwards 
into the halo. Figure [8] shows the constant (3 simulation 
DLEBc at the same age. The discussion of the differ- 
ences between the ExpCBc and ExpEBc applies also to 
simulations in the DL atmosphere. 

We performed an extra simulation (Tomis aka A) in 
a larg er volume to compare with Model A in ITomisakal 
(jl998|), with the same equatorial density. Other initial 
conditions were set as closely as possible to those spec- 
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Fig. 6. — Simulation DLH (hydrodynamic) at an age of 10 Myr. 
The density in a slice through the cube at the location of the 
source, (xi,X2,X3)=(0,0,0), in the xiX2 plane is shown with lin- 
ear grayscales from —0.2 (white) to 1.4 (black) cm~^. The unit of 
length on the axes is 100 pc. 



ified by iTomisakal (|1998D (see Tables [3] & i]) , but some 
small differences in the setup could not be avoided. Our 
simulations show the same characteristics in terms of the 
size and shape of the cavity, and the thickness of the sur- 
rounding shell. However we found differences of the order 
of 20% in the dimensions of the cavity after 10 Myr which 
appear to be associated with small differences in initial 
conditions. The fast magnetosonic shock had proceeded 
to almost twice the height of the shock in Model A of 
iTomisalal ()1998[ ). indicating a lower value of /3 in our 
simulation. The shape of the magnetosonic shock front 
is particularly sensitive to small differences in f3 because 
the Alfven speed at high altitudes is very sensitive to (3. 



4.4. Bubble Morphology as a Function of Time and 
Magnetic Field Strength 

Our simulations cover a factor 30 in magnetic field 
strength, from relatively weak (those labeled with 'a', 
f3 — 10) to strong (those labeled with 'd', /3 = 0.3). Here 
we discuss variation in the shape of a bubble with mag- 
netic field strength at a reference age of 10 Myr. The 
weak magnetic field limit for the simulations converges 
to the hydrodynamic case as should be expected. The 
differences between simulations with /3 = 10 and hydro- 
dynamic simulations are minor. We limit the discussion 
to the range /3 = 10 to 0.3. 

One of the most visible effects of the magnetic field is 
on the size of the cavity. A decrease in /3 by a factor of 30 
results in a decrease in the vertical size of the cavity by 
a factor ~1.7, independent of the density and magnetic 
field stratifications. The diameter of the cavity in the xi 
direction (along the field lines) increases by a factor of ^ 
1.4 as /? decreases. Expansion along the magnetic field is 
faster as the magnetic field is stronger. The diameter of 
the cavity along the X3 axis decreases by a factor ~1.6 
as (3 decreases from 10 to 0.3. The effect of a stronger 
magnetic field on the shape of the cavity is that the cavity 
becomes more elongated along the Galactic plane when 
observed along a direction perpendicular to the Galactic 
magnetic field. When observed along the direction of the 
magnetic field, the shape of the cavity varies less, but the 
size of the cavity becomes significantly smaller. 

The thickness of the shell in the directions perpendicu- 
lar to the magnetic field also increases dramatically with 
increasing magnetic field strength. However, this hap- 
pens at the expense of the amount of compression of gas 
in the shell, which becomes nearly invisible in the simu- 
lations with /3 = 0.3. In practice, the shell may not be 
detectable in real data confused with emission in the fore- 
ground and background, even for /3 = 1, if the magnetic 
field scale height is large (see Figure S]). The increasing 
thickness of the shell along the X3 axis is partly the re- 
sult of the smaller extent of the cavity, partly because of 
a faster expansion of the outer shock, approximately by 
equal amounts. In simulations with a constant magnetic 
field, the thickness of the shell at the top of the bubble 
increases mainly because of the fast upward expansion of 
the outer shock, and to a lesser amount by the smaller 
extent of the cavity. In the equipartition simulations, the 
increased thickness of the top of the shell is mainly be- 
cause of the smaller extent of the cavity. Along the mag- 
netic field lines (xi), the thickness of the shell decreases 
as the magnetic field strength increases. However, no 
strong increase in the density is found at the extremes of 
the shell along the xi axis. 

Figure[9]compares the time evolution of the semi- major 
axis of the contact discontinuity and the outer shock in 
the ii^a plane, for a sample of our simulations (ExpH, 
DLH, ExpCBc, & DLCBc), with that of the lGastor et all 
(|1975. ) model. The semi-major axis is defined as half the 
diameter along the xi axis at the level of the source. The 
expansion along the magnetic field lines at early times 
(see Figure |9K) close ly follows the hydrody namic solution 
in accordance with Ferrier e et al.l (|1991l ) . The contact 
discontinuity of the self-simila r solution was found to be 
at 86% of the shock radius by I Weaver et al.l (|1977D . The 
radius of the contact discontinuity in the analytic model 
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Fig. 7. — Simulation DLCBc at an age of 10 Myr. Panels show slices through the cube at the location of the source, (a::i,a;2,a;3)=(0,0,0), 
in three orthogonal planes. Grayscales show the gas density on a linear scale from —0.2 (white) to 1.4 (black) cm~^. The vector field 
depicts the projection of magnetic field vectors on the plane of the image. The unit of length on the axes is 100 pc. 
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. — Same as Figure [71 but for Simulation DLEBc at an age of 10 Myr. 
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indeed agrees well with the radii of the contact disconti- 
nuity in our simulations at early times (dashed curve in 
Figureini). ^t later times (FigurelSJs), the contact discon- 
tinuity in the ExpCBc and DLCBc simulations has pro- 
ceeded further than the analytic solution (dashed curve). 
On the other hand, the contact discontinuity in the ExpH 
and DLH simulations lags behind the analytic solution. 
We interpret this behavior as the result of the hydrody- 
namic simulations breaking out of the disk, reducing the 
pressure inside the cavity, whereas the MHD simulations 
remain confined to the disk. We find that the location of 
the outer shock in these simulations is nearly the same 
even after 15 Myr in agreement with results presented in 
lFerriereet"aII(|1991h . 

4.4.1. Axial Ratios 

The discussion in the previous subsection shows that 
the magnetic field is as important as the density dis- 
tribution of the medium in determining the shape of a 
supcrbubble. Observationally, the shape of the cavity is 
most important because structures with a large contrast 



in density are most easily detected in H I images of the 
Galaxy. 

Figure [10] shows the axial ratio of the low-density cav- 
ity in the Galactic plane (x^/xi) versus time for different 
resolutions (100'^,200"^ and 300"^) to verify convergence 
(see Section 21). From this, we confirm that convergence 
is established at 200^^ which we use for the analysis. Fig- 
ure [11] shows the relationship between the x^/xi axial 
ratio and (3 at an age of 5 Myr and 12.5 Myr. The simi- 
larity between the four panels shows that the axial ratio 
of a superbubble in the Galactic plane does not depend 
significantly on the vertical stratification of the density 
and the vertical stratification of the magnetic field (see 
also Figures [l] [5] [7] and [8]) . The important parameters 
are the age of the bubble and the strength of the mag- 
netic field in the equatorial plane. Even at the relatively 
young age of 5 Myr, a bubble in a medium with (3 = \ 
will be elongated along the magnetic field with axial ratio 
~0.8. After 12.5 Myr, the size of the bubble perpendic- 
ular to the magnetic field is only 60% of its size along 
the magnetic field. Variation of /3 by a factor of three 
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either way changes the predicted shape of the bubble in 
the Galactic plane from an axial ratio ~ 0.4 to ~ 0.8 at 
12.5 Myr. 

Figure [T^ shows the time evolution of the X2lx\ axial 
ratio of the low density cavity for ExpH (diamonds) , Ex- 
pCBc (stars), ExpEBc (filled squares), and the Kompa- 
neets model (dashed line). At early times, i ^ 5 Myr, the 
evolution of ExpH follows the Kompaneets model. The 
difference between the ExpH simulation and the Kom- 
paneets model at later times is the result of inertia of 
the swept-up interstellar medium. ExpCBc and ExpEBc 
evolve in a /3 = 1 medium which resists expansion along 
the xi axis. Contrary to ExpH, at early times, the axial 
ratios of ExpCBc and ExpEBc decrease, indicating that 
the expansion along the x\ axis is faster than along the 
xi axis. This trend continues for the ExpCBc simula- 
tion while it reverses for the ExpEBc simulation where 
the top of the bubble accelerates at later times. In the 
ExpEB simulations, confinement in the xi direction by 
the magnetic field decreases rapidly as the bubble grows 
beyond approximately one scale height for all values of 
/3. This allows the top of the bubble to accelaterate at 



Log(,S) 



Fig. 11. — Axial ratios (x-ijxx) of simulated bubbles as a func- 
tion of /3 at an age of 5 Myr (solid curve) and 12.5 Myr (dashed 
curve). In each panel dots represent the axial ratio measured in 
our simulations. The curves in all four panels converge to an ax- 
ial ratio of unity in the limit /3 — > oo (the hydrodynamic case). 
(A) Exponential atmosphere with constant magnetic field (Simu- 
lations ExpCBa-d), (B) Exponential atmosphere with constant /3 
(ExpEBa-d), (C) Dickey and Lockman atmosphere with constant 
magnetic field (DLCBa-d), (D) Dickey and Lockman atmosphere 
with constant /3 (DLEBa-d). The axial ratio at the beginning of 
the simulations is unity for each fi. 



later times, eventually leading to blow-out if the source is 
sufficiently strong (Tomisaka 1998). Breakout in our Ex- 
pEBc simulation occurs at t > 17.5 Myr. The constant 
magnetic field in all the ExpCB simulations confines the 
bubble in the X2 direction, preventing blow-out. 

Taking into account differences in density and luminos- 
ity of t he source, our axial ratios a gree fairly well with 
those of [ To misaka fc Ike uchil ()1986( ). while the results of 
iMac Lowet al.. (1989.') agree better with the Kompaneets 
solution. Differences between our results and those of 
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Fig. 12. — Axial ratio x^jxx as a function of time for a bubble in 
an exponential atmosphere with /3 = oo (diamonds; ExpH), /3 = 1 



with B ' 



ol/2 



(filled squares; ExpEBc), and /3 = 1 with constant 



B (stars; ExpCBc). Also shown is the axial ratio of the K60 model 
as a function of time (dashed curve). 
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Fig. 13. — Axial ratios X2/X1 as a function of time for a bubble 
in the Dickey and Lockman atmosphere with /3 = 00 (diamonds; 



run DLH), (3 = 1 with B ' 



»l/2 



/3 = 1 with constant B (stars; run DLCBc). 



(filled squared; run DLEBc), and 



iMac Low et al.l (119891 ) appear to be related to differences 
in source luminosity and the density distribution in the 
atmosphere. 

Figure [13] shows the time evolution of the X2/a;i axial 
ratio of the low density cavity for DLH (diamonds) , DL- 
CBc (stars), and DLEBc (filled squares). The evolution 
of the xijxi ratio of DLH shows the bubble accelerating 
in the X2 direction. The axial ratio at any time is larger 
than those measured for ExpH because the bubble is ex- 
panding in both the positive and negative xi directions. 
DLCBc and DLEBc evolve in a /3 = 1 magnetic field. 
The shape of the cavity in the x\X2 plane remains nearly 
circular for both magnetic field configurations, even at 
later times. A stronger magnetic field (/3 < 1) results in 
a cavity elongated along the xi-axis. As for the expo- 
nential atmosphere, the bubble in the DLEB simulations 
eventually breaks out, whereas the bubbles in the DLCB 
simulations remain confined by the magnetic field. 

4.4.2. Deriving Scale Height and Age From Analytic Models 

The elongation of a bubble along the magnetic field has 
important consequences for age estimates of the bubble 



from observations, which usually assume axial symmetry 
or even spherical symmetry. The axial ratio x^jxi can be 
observed for superbubbles in face-on galaxies, but not for 
Galactic superbubbles which are observed from within 
the Galactic plane (2:1X3 plane). Within the Galactic 
plane, the available observations are the radius of the 
bubble in the plane of the sky and the expansion veloc- 
ity perpendicular along the line of sight. The age of the 
bubble derived from observations is to\,s = OiRobs/vexp-^ 
assuming an expansion law of the form R '^ t"'. The 
elongated shape of the bubble along the direction of the 
magnetic field implies that, on average, the expansion ve- 
locity is larger along the xi-axis than along the X3 axis. 
An observer who looks at the bubble along the x\ axis 
(looking along the magnetic field lines) will see the small 
dimension of the bubble along the x^ axis but measure the 
larger expansion velocity along the xi axis. An observer 
who looks at the bubble along the x^ axis (looking per- 
pendicular to the magnetic field lines) will see the large 
dimension of the bubble along the xi axis but measure 
the smaller expansion velocity along the x^ axis. The 
two observers will derive significantly different ages for 
the bubble from their observations. If the radius of the 
bubble along the Xi axis is i?o and the axial ratio x^/xi 
is <7, the age derived by observer A (looking along the 
field lines) is 

t^ = a-^^, (22) 



w, 



exp^A 



and the age derived by observer B (looking perpendicular 
to the field lines) is 



ts = Oi- 



Rq 

Vexp,B 



(23) 



On average, Vexp,B — qVexp,A, therefore the ratio of the 
ages derived by the two observers is 



tA 
tB 



(24) 



For a moderate x^/xi axial ratio of 0.8 (Figure [TT]), the 
ratio tAJtB = 0.64. Age estimates for observer A looking 
along the Xi axis, and observer B looking along the 2:3 
axis were derived from our simulations. The ratio of 
these age estimates is even somewhat smaller than q^, 
especially for lower values of (3, because the expansion of 
the shell actually stalls in the X3 direction after a finite 
amount of time. This illustrates the importance of MHD 
effects on the derivation of the age of a superbubble. 



TABLE 5 
Fitting the Kompaneets Model to Simulations 

IN the X3X2 PLANE AT t = 7.3 MYR 
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The three-dimensional structure of a magnetized bub- 
ble not only affects age estimates, but also estimates of 
the scale height of the surrou nding medium from the ob- 
served shape of the bubble. iBasu et ahl ()1999f ) outline 
a method for fitting the Kompaneets model to observa- 
tions to obtain the scale height of the ambient medium 
and the age of a bubble. Our simulations can be used 
to assess systematic errors introduced by applying an 
axially symmetric hydrodynamic model to a magnetized 
super bubble. Inspection of Figures IH [H [71 and [5] shows 
that in our simulations the Kompaneets model can only 
fit the shape of the cavity in an exponential atmosphere 
looking approximately along the magnetic field lines, or 
the circular shape of the cavity in a DL atmosphere look- 
ing approximately perpendicular to the field lines. In the 
latter case, fitting the Kompaneets model to the circular 
shape of the cavity would imply a very large scale height. 
Here we focus on fitting the Kompaneets model to the 
ExpH, ExpCBc, and ExpEBc simulations, in the 2:3X2 
plane. 

To determine a range of Kompaneets solutions that 
fit these bubbles, an initial estimate for y (Equation [4]) 
was derived from the xzjxi axial ratio of the simulation. 
Starting with this estimate, overlays of the Kompaneets 
solution were visually fitted to the simulations to obtain 
a range of acceptable y values, with the additional con- 
straint that the source of the Kompaneets model must 
coincide with the source in the simulations to within 2 
pixels, i.e. the radius of the source in the simulations. 

Fits of the Kompaneets model have also be used to 
estimate the age of a bubble, if the mechanica l luminosity 
of the source and ambient density are known. iBasu et al.l 
(|1999f ) showed that at the level of the source, the radius 
of the Kompaneets bubble is nearl y equal to the radiu s 
of a spherical bubble described by I Cast or et al.l (|1975D . 
From Equation ([3]), at z = 0, we find 



R[z = 0, y) = 2iJ arccos 1 — 



y 



(25) 



As y and H are determined from the fit. Equation (j25p 
can be substituted into Equation ^ with Lj = 3 x 10'^^ 
erg s~^ and p = 1.67 x lO^^"' g cm^'^ to find the age of 
the bubble. 

Table [5] shows results of visual fits of the Kompaneets 
model to the ExpH, ExpCBc, and ExpEBc simulations at 
an age of 7.3 Myr, looking along the xi axis. Scale height 
and age are listed as a function of y for those models that 
provided an acceptable fit. For all three simulations, the 
fitted scale height is smaller for larger values of y, sim- 
ply indicating a range of acceptable axial ratios that fit 
the shape of the bubble. Simulations ExpH and ExpCBc 
are fitted by the same range of y while ExpEBc is better 
fitted by higher values of y. This reflects the more elon- 
gated shape of the cavity in ExpEBc seen in Figure [5l 
For ExpH, we find that the scale height is close to 100 
pc. Although ExpH and ExpCBc are fitted by the same 
range of y, the scale height derived from fitting ExpCBc 
is smaller because the bubble is more confined. This con- 
finement of the bubble by the magnetic field while main- 
taining an elongated shape of the cavity is the root cause 
of the 30% to 50% lower scale heights resulting from fits 
of the Kompaneets model to a magnetized bubble. Ages 
derived from fits to the magnetic simulations are consis- 
tently a factor ~ 2 smaller than the actual age of the sim- 



ulation. We note in passing that iKoo fc McKed (|1990[ ) 
found that the Kompaneets model blows out a factor ~ 2 
earlier than their numerical hydrodynamic solutions. 

As discussed earlier, the bottom of the bubble in the 
K60 model docs not penetrate deeper than 1.4 exponen- 
tial scale heights below the level of the source. If a bubble 
would penetrate deeper into the atmosphere, an observer 
could conclude that the Kompaneets model is not appli- 
cable. However, none of our simulations penetrate as 
deep as the limiting value of 1.4 scale heights below the 
level of the source. In practise, this would not consitute 
a conclusive test whether or not the Kompaneets model 
can be applied. 

5. APPLICATION TO THE W4 SUPERBUBBLE 

As a specific ap plication of our si mulations we investi- 
gated the result of lBasu et al.1 ()1999l ) that the scale height 
of the interstellar medium around the W4 superbubble 
(jNormandeau et anil996[ ) is as small as 25 pc, approx- 
imately a fact or 4 smaller than the canonical value of 
100 pc. iKomlienovic et a ll |l999) investigated this issue 
with two-dimensional MHD simulations, and suggested 
that the elongated shape of the bubble could only be gen- 
erated by a magnetic field perpendicular to the Galactic 
plane (see also West et al. 2007). Our simulations allow 
us to consider the elongated shape of the bubble seen 
by an observer looking along the magnetic field lines. 
This perspective is not possible in two-dimensional sim- 
ulations. This section is organized as follows. First we 
review some basic properties of the W4 region, because 
we found that the ambient density is a factor ~ 2.5 lower 
than adopted by previous authors. Second, we proceed 
w ith our analysis of the sca le height near W4. 

iNormandeau et ahl ()1996l ) showed that an expanding 
shell in the Perseus arm associated with the Cas 0B6 
association (Braunsfurth 1983) is a conical cavity with 
no apparent upper boundary in HI. This structure was 
interpreted as a "chimney" , a large stellar wind bubble 
which has broken out of the Galactic disk and expands 
into the Galactic halo. Evidence for a vertical fiow inside 
the cavity comes from el ongated H I structures ass ociated 
with a molecular cloud (INormandeau et al]|1996f), and a 
comet-shaped mo lecular cloud (jNormandeau et al.lll996t 
iHever et al.lll996[ ). The lower part of the bubble wall is 
brigh t in Ha and is also s een in radio continuum images 
(Nor mandeau et al]|1996D . The bright H 11 region at th e 
bottom of the cavity is known as W4 (|Westerhoutlll958t) . 
The bubble appears to be confined at the bottom by a 
dense cloud. 

The probable source of the bubble is the star cluster 
OCL 352 (IC 1805), which is located near the bottom 
of the chimney. INormandeau et al.l (jl996( ) list 9 O stars 
in the cluster, with the earliest spectral type 041 for the 
star BD-h60504. The presence of the 04 star in OCL 352 
makes it likely that no supernova explosion has yet oc- 
cured in the cluster, so that the bubble is a true stel- 
lar wind bubble. Most recent age determinatio ns for 
the cluster are in the range 1 to 3 Myr (.Hillwig et al.1 
l2006HRauw fc De Beckeill2004l:lMassey et al.lll995D. but 



ages up to 5 Myr have been suggested (|Kharchenko et aLJ 
|2005[) . The kinetic energy released in the stellar wind 
of the stars in the cluste r is Lmech = 3 x 10^^ ergs~^ 
(INormandeau et aLlll996t ). most of which is emitted by 
BD-F60504 and two 05 stars. 
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iDennison et al.l (|1997f ) reported an elongated Ha shell 
surrounding the HI cavity that appears closed 6° (230 
pc) above the star clust er, outside the field of vi ew of 
the original H I images of lNormandeau et al.1 (|1996D . Ex- 
tended HI images that cover the top of the Ha shell 
do n ot show a cap in the neutral gas (jNormandeaul 
I2OOOI) . The Ha shell was interpreted as the inner wall 
of the superbubble io nized by the stars in OCL 352. 
IDennison et al.1 (|1997[ ) found that the age of the super- 
bubble exceeds the maximum possible age of OCL 352 
by a factor ~ 3, using the mechanical luminosity of the 
star cluster {L„iech = 3 x 10'^^ ergs~^) and the den- 
sity of the ambien t med ium (no = 5 cm~'^) given by 
iNormandeau et al] (|1996D . The large dynamical age of 
the bubble was taken as an indication that the star clus- 
ter OCL 352 by itself may not b e the only source of the 
superbubble. iBasu et aLl (|1999D found that the age of 
the bubble is consistent with the age of the star clus- 
ter by applying the K60 model to the W4 superbubble, 
assuming a scale hei ght of only 25 pc for the ambient 
medium. However, iDennison et al.l (1997) noted that a 
lower ambient density could decrease the dynamical age 
derived for the W4 superbubble. 

The original H I number density was derived by 
INormandeau et aLl ()199 6l from an HI column density 
map. This map was obtained by integration of the 
HI brightness temperature over a velocity range of 
14.8 km s~^. The ambient density was then calculated 
by dividing the HI column density adjacent to the cavity 
by the diameter of the bubble, assuming that all the HI 
emission is local to the W4 superbubble. This assump- 
tion leads to a significant over-estimate of the ambient 
density, because H I emission in the velocity range of the 
bubble originates from a much longer section of the line 
of sight. We re-analyzed the HI data of the W4 region, 
now publicly available as part of the Canadian Galactic 
Plane Survey (JTavlor et al.ll2003f ). to re-determine the 
ambient density. 

The HI column density inside the cavity of the W4 
superbubble, avoiding the HI features known to be inside 
the cavity (Normandcau et al. 1996), is 1.1 x 10^^ cm~-^. 
The medium inside the superbubble is expected to be a 
highly ionized tenuous plasma at a temperature of a few 
million Kelvin. The ionized bubble wall seen on both 
sides of the cavity in radio continuum and Ha emission 
shows that the interior of th e bubble is optically th in to 
Lyman continu um photons. IDennison et al.l (|1997f ) and 
IBasu et "aLl ()199 9) assume this to be true in their analysis 
of the ionization of the bubble wall by the star cluster. 
Therefore, HI observed inside the cavity must be located 
in the foreground and the background. Some of this HI 
emission could in principle belong to the bubble wall if 
the expansion of the superbubble had stalled. There is 
no reason to believe that the expansion of the bubble has 
stalled because the source is still active. 

Although the evidence tha t the H I steamers identified 
bv INormandeau et al.1 ([1996) are located inside the cav- 
ity of the W4 bubble is compelling, the HI emission that 
seems to partly fill the cavity everywhere is likely unre- 
lated gas, that was inclu ded in the density estimate of 
INormandeau et aLl ()1996t ) . If the H I column density in- 
side the cavity is subtracted, the density of the ambient 
medium is reduced to ~ 2 cm~^ assuming the same line- 
of-sight dimension of the bubble. A better estimate of 



the density of the ambient medium is obtained by divid- 
ing the mass of the shell of swept-up gas by the volume 
of the cavity. This is impractical for the W4 bubble be- 
cause there is no clear limb-bri ghtened H I shell visible 
in the column density map. IBasu et al.1 (|1999[) found 
a higher density of 10 cm"'^ from their analysis of the 
shape of the ionization front, which depends on the de- 
tails of their model. Although the lower value of the am- 
bient density derived from the HI data may not resolve 
the age problem altog ether, we follow the suggestion by 
IDennison et al.1 (|199 7t) that the location of the cluster is 
strong evidence for its identification as the source of the 
superbubble. 

We no w return to the is su e of the smal l scale height 
found bv IBasu et al.l (|1999f) . iNormandeaiil (|2000) found 
no evidence for such a small scale height in HI images 
of the region, but, as noted before, the HI images are 
contaminated to some extent by emission in the fore- 
ground and background. The presence of molecular gas 
with a small scale hei ght may have contributed to the 
result obtained by Basu et al.l ([1999) , but the molecu- 
lar gas traced by the CO fine appears too fragmented to 
confine the shape of the bubble. The scale height found 
by fitting the Kompaneets model to the observed shape 
of the bubble may also be artificially small by a factor 
2 (see Ji4.4.2p for an observer looking along the magnetic 
field parallel to the Galactic plane. This factor is not 
enough to explain the very small scale height derived by 
IBasu et aLJ ()1999[ ). Also, the age derived bv IBasu et al.1 
|199sD would be too small by a factor 2, further increas- 
ing the discrepancy between the dynamical age of the 
bubble and the age of the star cluster. None of our sim- 
ulations evolve into a cavity that is as narrow as the W4 
super bubble. 

6. FARADAY ROTATION BY A MAGNETIZED BUBBLE 

Polarimetric observations of the Galact ic plane 
at arcminute reso l ution (jTavlor et al.l 120031 : 
iMcClure-Griffiths et al.1 l2005f ) allow measurements 
of Faraday rotation of diffuse Galactic synchrotron 
radiation and background compact extragalactic sources 
to probe the magneto-ionic interstellar medium down 
to parsec scales. The plane of polarization of linearly 
polarized synchrotron emission rotates by an angle A0 
(radians) as the radiation passes through a magnetized 
plasma according to /S.9 = X^RM, with A the wavelength 
of the radiation in meters. The line-of-sight integral 



RM = 0.812 / UeB ■ dl 



(26) 



is the rotation measure (radians m^^), with Ug^ the elec- 
tron density (cm~^), and B • dl the magnetic field vec- 
tor B (/iC) projected onto the line of sight (in parsec). 
The rotation measure can be determined observationally 
by observing A9 at different, but closely spaced wave- 
lengths. Structure along the line of sight may create a 
measurable rotation of the plane of polarization, but it 
may not have a sufficiently large emission measure for 
any thermal emission to be detected. 

Equation (PS|) cannot be inverted to obtain a unique 
solution to the magnetic field structure, even if indepen- 
dent information on the electron density Ue along the 
line of sight is available. Instead, the data can be in- 
terpreted by line-of-sight integration of a model distri- 
bution of the density and magnetic field, which depends 
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Fig. 14. — Images of rotation measure of the simulation DLCBc at age 10 Myr (see also Figure (T)! along a direction perpendicular to 
the undisturbed Galactic magnetic field (left) and along the direction of the magnetic field (right). Grayscales are linear from —20 to +20 
radians m~^ (left) and from 130 to 3252 radians m~^ (right). Wrapping of the magnetic field lines around the low-density cavity creates 
significantly different Faraday rotation signatures for the same bubble seen from different directions. If the line of sight is perpendicular to 
the Galactic magnetic field, the shell of swept-up interstellar medium is almost invisible in the rotation measure map. If the line of sight 
is along the Galactic magnetic field, the largest rotation measures are associated with the (thick) shell. The unit of length on the axes is 
100 pc. 

general characteristics of the Faraday rotation of all bub- 
bles in our simulations. As the magnetic field is pushed 
aside, the field lines wrap around the expanding cavity. 
The largest amplification and the largest change of direc- 
tion of the magnetic field occur just outside the cavity, 
in the equatorial plane (for the DL density distribution). 
The density is also enhanced, compared with the sur- 
rounding undisturbed interstellar medium, but by differ- 
ent amounts depending on location in the shell. Sim- 
ulations with a small scale height of the Galactic mag- 
netic field (B ~ p^^^) show this compression also at the 
top of the bubble, whereas simulations with a large scale 
height of the magnetic field (constant B) do not. The 
largest Faraday rotation therefore occurs in a relatively 
thin region around the cavity, being strongest close to the 
Galactic plane (see Figure [H)l. The Faraday rotation in- 
side the cavity is much smaller by comparison, because of 
the low density and the low chaotic magnetic field there. 
If the line of sight is perpendicular to the direction of 
the undisturbed magnetic field, the near and the far side 
of the bubble wall have rotation measures up to ~100 
radians m^^, with opposite signs. These mostly cancel 
each other, but asymmetries between the front side and 
the back side of the shell create structure in the rotation 
measure distribution. Asymmetries arise from the insta- 

The 



critically on the assumed distributions of the electron 
density and magnetic field. If an ad hoc model for den- 
sity and magnetic field is adopted to interpret the ob- 
servations, careful consideration should be given to ba- 
sic questions regarding self-consistency, uniqueness and 
dynamical properties of the proposed model. Our three- 
dimensional MHD simulations of magnetized superbub- 
bles constitute a significant step forward in the modeling 
of Faraday rotation by Galactic superbubbles. In this 
paper, we only consider the bubble as a Faraday screen 
that does not emit synchrotron emission. This situation 
applies in particular to compact polarized sources in the 
background. 

The simulations provide a self-consistent set of mod- 
els that can be used to calculate the rotation measure 
imposed on background emission by a magnetized super- 
bubble. The non-spherical evolution of the superbubble 
described in Section |4] may have very significant effects 
on the predicted rotation measure. The effect of nearby 
superbubbles on Far aday rotation of background sources 
was first discussed bv lValleg ()1993 ^. In this paper we dis- 
cuss only the Faraday rotation signature of distant bub- 
bles, for which we can assume that all lines of sight are 
parallel. We assume here that the bubble is completely 
ionized. This allows a first exploration of the general ef- 
fects of the topology of the density and magnetic field, 
without introducing more parameters. 

Figure [14] shows the rotation measure for a line of sight 
along the X3 axis (looking perpendicular to the Galactic 
magnetic field) and the Xi axis (looking along the Galac- 
tic magnetic field) of the simulation DLCBc (Table [3]) at 
an age of 10 Myr (Figure[7]). The linc-of-sight integration 
of Equation (|26p was done through the entire simulated 
volume as shown in Figure [T] It includes some undis- 
turbed medium in the foreground and the background of 
the bubble. 

Comparing the rotation measure maps with the den- 
sity and magnetic field structure (Figure [7]) reveals some 



bilities in the shell mentioned in Sections l4.2.2l fc l4r41 
strongest signal is expected from asymmetries that occur 
in locations where a perfectly symmetric bubble would 
have bent the Galactic magnetic field sufficiently to cre- 
ate a significant line of sight component of the magnetic 
field. As the magnetic field is tightly wrapped around 
the cavity, the highest values of the rotation measure 
resulting from front-to-back asymmetries are therefore 
expected for lines of sight that intersect the cavity, away 
from the vertical axis that intersects the source. Our 
three-dimensional simulations contain such asymmetries, 
resulting in the structures seen in Figure [14] with ampli- 
tudes of the order of 20 radians m^^. 
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Fig. 15. — Simulation ExpCB with cooling at an age of 10 Myr for comparison with Figure|4] Magnetic field initially oriented along xi 
axis, (3 = 1. Panels show slices through the cube at the location of the source, (xi,X2,a;3)=(0,0,0), in three orthogonal planes. Grayscales 
show the gas density on a logarithmic scale from 10~^ (white) to 10^'^ (black) cm~^. The vector field depicts the projection of magnetic 
field vectors on the plane of the image. The unit of length on the axes is 100 pc. 
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Fig. 16. — Simulation DLCB with cooling at an ago of 10 Myr for comparison with Figure[7] Magnetic field initially oriented along xi 
axis, (3 = 1. Panels show slices through the cube at the location of the source, (a;i,a;2,a;3)=(0,0,0), in three orthogonal planes. Grayscales 
show the gas density on a logarithmic scale from 10~^ (white) to lO'^'^ (black) cm~^. The vector field depicts the projection of magnetic 
field vectors on the plane of the image. The unit of length on the axes is 100 pc. 



If the line of sight is paraUel to the direction of the 
undisturbed magnetic field, rotation measures are much 
higher everywhere than in the previous case. The front 
and back side of the shell reinforce each other, and high 
rotation measures are expected for lines of sight that in- 
tersect the shell. The top of the bubble increases the ro- 
tation measure a few hundred parsecs above the Galactic 
plane by ~ 20% compared with the undisturbed atmo- 
sphere in Figure [T31 For the simulation DLEBc (small 
magnetic scale height), the increase is ~ 30% compared 
with the undisturbed atmosphere. Within 100 pc from 
the Galactic plane, the rotation measure varies with lo- 
cation because of the bubble, but the mean of the rota- 
tion measure across the bubble is almost the same as the 
mean rotation measure of the undisturbed atmosphere. 
Here, the effect of the bubble is to create a significant 
variation in rotation measure for different lines of sight, 
with little effect on the mean rotation measure taken over 



a large area. Super bubble shells may therefore enhance 
the rotation measure in the disk-halo interface and affect 
estimates of the scale height of the Galactic magnetic 
fi eld from rotation measu re analysis. 

IVallee fc Bignelll ()1983D reported a large magnetized 
bubble associated with the Gum nebula, based on a 
number of high ro tation measures towa rds extragalactic 
sources. Recently, ' Stil fc Tavloil {2Wf} found this bub- 
ble in an all-sky image of depolarization of extragalactic 
sources, that was associated with excess rotation mea- 
sure. This bubble stands out in the data because of the 
high rotation measure in the shell. The line of sight in 
this case is nearly parallel to the magnetic field in the 
local spiral arm. 

7. THE EFFECTS OF COOLING 

So far we have discussed adiabatic simulations of super 
bubbles as an approximation that energy losses are com- 
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pensated by an internal or an external radiation field. 
In this section we briefly discuss simulations with cool- 
ing in the two main density stratifications discussed in 
this paper: the exponential atmosphere and the Dickey 
& Lockman atmosphere. The cooling represents an extra 
term in the energy equation (10) and is implemented as 
an explicit source term in the code. As with all physical 
processes, there is a maximum allowed time step for nu- 
merical stability associated with the cooling; in our case 
it is reset to 10% of the radiative cooling time when it is 
below the normal hydrodyn amic CFL condi ti on. 

The cooling funct io n by iRavmond et al.1 ()1976l ) and 
Dalgarno fc: McCravl (|1972l ). using the abundances of 



AllenI (|1973l ). was^implemented into zeus-mp in the 



form of analytic expressions adopted by 'Tomisak a et al.l 
([198 li) . Although some abundances were revised down- 
wards later, the difference between the cooling rates 
adopted in ou r simulations and t h e sola r abundance cool- 
ing curve of iGnat &: Sternberg ()2007f ) is smaller than 
other sources of uncertainty introduced by finite reso- 
lution, and ionization by the central star cluster. If an 
accuracy better than a factor ^ 3 is required in the cool- 
ing rate, the location of the superbubble in the Galaxy 
also becomes important because of variati on in the abun- 
dances with location in the Galaxy (e.g. iRudolph et al.l 
[2006). 

The radiatively cooled simulations were done in atmo- 
spheres maintained at a temperature of 8000 K, with the 
same resolution as the adiabatic simulations discussed 
in Section IH and with the same mechanical luminosity 
of the source. In these simulations, cooling begins at a 
time t = A Myr, soon after the first supernova explo- 
sion, at the beginni n g of t he radiative phase following 
iMacCrav fc Kafatod l)l987D . The effects of cooling on 
the shape of the cavity are evaluated at i = 10 Myr. 

7.1. Results with cooling 

Figure [15] and Figure [12] show simulations of a bub- 
ble with cooling in an exponential atmosphere and in a 
Dickey & Lockman atmosphere at an age og 10 Myr. 
These figures can be compared with the adiabatic sim- 
ulations shown in Figure [5] and Figure [7] Qualitatively 
the features in these simulations are very similar. The 
cavity is elongated along the direction of the magnetic 
field. The shell is thicker in directions perpendicular to 
the magnetic field in the undisturbed atmosphere {x2 
and xs), while significant compression of the swept-up 
material in the direction of the magnetic field leads to a 
thinner, denser shell. Closer inspection shows that the 
cavities of both cooled simulations are more elongated in 
the direction of the magnetic field [xi). The shell also 
seems to display stronger instabilities, leading to irregu- 
larities in the shape of the cavity, and increased scatter 
in the derived axial ratios of the cavity. As expected, 
the simulations with cooling lead to shells that are thin- 
ner, with a higher density and stronger magnetic field, 
than the adiabatic simulations. The effect of resolution 
on the cooled regions is to smear out small-scale high- 
density structure wit h the shortest cooling time and the 
lowest temperature. iDe Avillez fc Breitschwerd^ ([2004f ) 
found convergence in the maximum density and lowest 
temperature in their simulations for a grid size ^1.1 pc. 
The density and temperature inside the swept-up shell in 
our simulations are probably also affected by resolution. 
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Fig. 17. — Axial ratio x^jxx (A) and x^/xi (B) as a function 
of the age of the bubble for cooled simulations in an exponential 
atmosphere (stars) and a Dickey & Lockm an a tmo sph ere (squares). 
Panel (A) may be compared with Figures Il2l and [131 and panel B 
may be compared with Figure [TOl 

but this paper focuses on the s hape of the bubble, not 
the co nditions inside the shell. iGarcia-Segura fc Loped 
([2000f ) found that the shape of planetary nebulae in their 
simulations was not affected by resolution. 

Figure [T7] shows the evolution of the axial ratios as a 
function of time in these simulations. We find that the 
axial ratio xijxx of a bubble in the Dickey fc Lockman 
atmosphere is most affected by cooling (compare Fig- 
ure [TTk with Figures [13] and [T3|) . In the adiabatic sim- 
ulation, the cavity remains nearly circular to an age of 
10 Myr, while in the cooled simulation the cavity is more 
elongated in the direction of the magnetic field. Note the 
change in evolution of the x^jxx axial ratio just after the 
onset of cooling. The effect of cooling on the axial ra- 
tio X2I x\ compared with the adiabatic case is less strong 
for the exponential atmosphere, but here too the cavity 
becomes more elongated in the direction of the magnetic 
field. A stronger elongation is also seen in the x^/xi ax- 
ial ratio (compare Figure 117b with Figure llOp , but the 
effect is smaller than on the axial ratio xijxx. 

The observed change in axial ratios in the presence of 
cooling is most likely the result of increased magnetic 
tension exerted by the stronger magnetic field in the 
dense shell that surrounds the cavity. This restricts the 
expansion perpendicular to the field (x2 and ^3) more 
in the cooled simulations than in the adiabatic simula- 
tions. Coupling between he magnetic field and the neu- 
tral shell is maintained through a small ionized fraction 
in the shell, and ion-neutral collisions. A complete pa- 
rameter study of simulations with cooling is beyond the 
scope of this paper. 

At least some superbubbles have a substantial fraction 
of the mass of the shell ionized by the stars inside. In 
simulations with cooling, the cooling time of the shell 
is very short compared to the age of the bubble, and a 
dense neutral shell forms. Ideally, radiative transport of 
the ionizing flux of the star cluster and any unrelated 
massive hot stars in the simulation volume should be 
included in the simulations. Solving the radiative trans- 
port every time step in the simulations is a computational 
challenge. The mean free path of an ionizing photon is 
substantially smaller than the resolution of the Simula- 
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tions. The ionizing flux decreases rapidly with time as 
the most massive stars explode as a supernova after ~ 3 
Myr, and the shell will eventually cool and become neu- 
tral. The implicit assumption in this paper by applying 
adiabatic evolution to this problem is that heating by 
photo-ionization is balanced by cooling. While this is 
a simplification, the observed ionization of superbubble 
walls suggests that heating by ionization of the shell is 
a non-negligible term in the energy equation of the shell 
that is not included in simuations with cooling only. 

8. CONCLUSIONS 

We present three-dimensional MHD simulations of a 
superbubble evolving in a magnetized medium. In these 
simulations, we assume two different atmospheric mod- 
els, exponential and Dickey & Lockman (1990), and two 
different magnetic field configurations, constant mag- 
netic field and constant /3, for varying values of the mag- 
netic field strength. With these simulations we aim to 
study the importance of MHD effects on the interpreta- 
tion of observed super bubbles. 

As noted before by iTomisakal (jlQQSD . a superbubble 
in a magnetized medium becomes significantly elongated 
along the magnetic field. We find that the axial ratio of 
the bubble in the Galactic plane is ^ 0.8 after 5 Myr and 
~ 0.6 after 12.5 Myr, depending on the age of the bub- 
ble and the strength of the magnetic field but not on the 
vertical structure of either the magnetic field or the gas. 
The elongated shape of the bubble in the Galactic plane 
may lead to significant errors in the age determination of 
magnetized superbubbles from observations using sym- 
metric hydrodynamic models. The derived age depends 
on the location of the observer. The ratio of the smallest 
age to the largest age derived by observers assuming ax- 
ial or spherical symmetry is found to be proportional to 
the square of the axial ratio of the bubble in the Galactic 
plane, which creates a discrepancy in the age of up to a 
factor ~4. 

We have analyzed systematic errors in the age of the 
bubble and the scale height of the ambient medium intro- 
duced by fitting the Kompancets model to a magnetized 
superbubble looking along the magnetic field lines. The 



scale height may be underestimated by 30% to 50% and 
the age by 50%. In particular, we investigated the curi- 
ously small scale height of th e interstella r medi um near 
the W4 superbubble found bv lBasu et al.l (|1999D . We re- 
analyzed HI data from the Canadian Galactic Plane Sur- 
vey and found that the density of the ambient medium 
uh ~ 2cm~^ which is smaller than previously thought. 
This lower density helps to diminish the apparent age 
discrepancy between the W4 bubble and the star cluster 
OCL 352. However, our analysis of the systematic errors 
introduced by fitting the Kompaneets model shows that 
they are not big enough to exp lain the scale height of 25 
pc found bv lBasu et al.l ()1999( ) for the W4 region. 

We use the MHD simulations to predict the rotation 
measure distribution of superbubbles based on three- 
dimensional MHD simulations, and emphasize the im- 
portance of such simulations to make these predictions. 
As expected, the appearance of a magnetized superbub- 
ble depends on the perspective of the observer. If the 
observer looks along the magnetic field lines, the largest 
rotation measures are seen at the intersection of the shell 
with the Galactic plane. The rotation measure is in- 
creased at larger distances from the Galactic plane. If 
an observer in the Galactic plane looks perpendicular 
to the magnetic field lines, the rotation measures are 
much smaller, but most importantly most of the struc- 
ture in rotation measure appears in projection on the 
low-density cavity, and not on the shell surrounding it. 

The simulations and analysis presented in this paper 
highlight the importance of three-dimensional MHD sim- 
ulations of superbubbles evolving in the Galactic mag- 
netic field to the interpretation of new high-resolution 
images of the Galactic plane at radio wavelengths from 
the International Galactic Plane Survey. 
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APPENDIX 
THE KOMPANEETS SOLUTION AT EARLY TIMES 

We show that the Kompaneets solution at early times can be approximated as an expanding sphere with a center 
that moves upward from the location of the source. This analytic derivation explains quantitatively the behavior of 
the Kompaneets solution in Figure [3l 

The center of the bubble along the vertical axis follows from the expressions for the top and bottom of the bubble 
(Equation [6]) 

,, = £l±ii = -Hln[l-(y/2f]. (Al) 

The half-diameter of the Kompaneets solution in the vertical direction is 

Zl - Z2 



Rz 



h + y/2^ 



(A2) 



T he maximum radius Rh of the Ko mpaneet s solution in a direction perpendicular to the density gradient was given 
bv lBisnovatvi-Kogan fc Silichl (|1995f ) and .Basu et al.l (|1999D . 



Rh = 2i?arcsin(y/2). 
The Taylor expansion of z^, Rz, and R^ in x = y/2 to third order in x is 



(A3) 
(A4) 



3D Simulations of Superbubbles 



21 



where 0{x ) is the remainder that contains terms of order a;* or higher. Similarly, we have 



Rz 



2Hx+-Hx^ + 0{x^), 



(A5) 



and 



so the difference R, — Ri, is of order x^ 



Rh 



2Hx 



l""' 



1 



R^ — Rh — ~Hx 



O(x^), 



0(x5). 



(A6) 



(A7) 



We see that for small x, the shift in the centre is Hx^ , which is of second order in x. The difference Rz — Rh is of third 
order in x. The Kompaneets solution for small y can therefore be approximated by a spherical bubble that rises in 
the atmosphere. The center according to equation lAll is shown in Figure [31 The offset is clear, where the Kompaneets 
model is still visually a circle. Although both the numerical simulation and the Kompaneets model have an axial ratio 
that is near unity, the difference between the Kompaneets model and the numerical simulation is in the fact that the 
former is not centered on the source. 
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